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Abstract — We propose a fast approximate algorithm for large graph matching. A new projected fixed-point method is defined and a 
new doubly stochastic projection is adopted to derive the algorithm. Previous graph matching algorithms suffer from high computational 
complexity and therefore do not have good scalability with respect to graph size. For matching two weighted graphs of n nodes, our 
algorithm has time complexity only 0{n^) per iteration and space complexity 0{n'^). In addition to its scalability, our algorithm is 
easy to implement, robust, and able to match undirected weighted attributed graphs of different sizes. While the convergence rate 
of previous iterative graph matching algorithms is unknown, our algorithm is theoretically guaranteed to converge at a linear rate. 
Extensive experiments on large synthetic and real graphs (more than 1 ,000 nodes) were conducted to evaluate the performance of 
various algorithms. Results show that in most cases our proposed algorithm achieves better performance than previous state-of-the-art 
algorithms in terms of both speed and accuracy in large graph matching. In particular, with high accuracy, our algorithm takes only a 
few seconds (in a PC) to match two graphs of 1 ,000 nodes. 

Index Terms — Graph matching, projected fixed-point, large graph algorithm, feature correspondence, point matching 
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1 Introduction 

GRaph matching, aiming to find the optimal cor- 
respondences between the nodes of two graphs, 
is an important and active topic in computer vision 
and pattern recognition [IJ, [2J. It has been extensively 
applied in various fields including optical character 
recognition |i3J, ||4J, object recognition [5], [6], shape 
matching O, ||7|, O, face recognition |9|, feature corre- 
spondence [lOj, point matching image retrieval |[T2l , 
and video indexing [13]. 

Graph matching is in general an NP-hard discrete opti- 
mization problem. Exact graph isomorphism algorithms 
include Ullman's method [14J, Nauty [15] and VF2 [16], 
all of which have exponential time complexity in worst 
cases. To match two graphs within a reasonable time, 
one has to look for approximate solutions. Moreover, due 
to noise and variability in real world graphs, the usage 
of exact graph matching algorithms is very limited. 
The focus of this paper is the design of a approximate 
algorithm for efficiently matching general large graphs 
(e.g. graphs of 1,000 nodes) in computer vision and 
patter recognition. 

One approach of approximate graph matching algo- 
rithms is based tree search ITTI , [18 1. Its basic idea is tree 
search with backtracking while using heurstics to prune 
unfrutful paths. Another approach of graph matching 
algorithms is based on continuous relaxation of the dis- 
crete problem while using cotinuous optimization tech- 
niques or heuristics to optimize a matching score. Classic 



• Yao Lu, Kaizhu Huang and Cheng-Lin Liu are with National Laboratory 
of Pattern Recognition, Institute of Automation, Chinese Academy of 
Sciences 

E-mail: yaolubrain@gmail.com, kzhuang@nlpr.ia.ac.cn, liucl@nlpr.ia.ac.cn 



work include Relaxation Labeling (191, EQ), ED, Ea 
and Graduated Assignment |[23l . Generally, continuous 
relaxation based algorithms have lower computational 
costs than heurstic search based ones [1]. 

In this paper, we propose a novel fast graph matching 
algorithm called Fast Projected Fixed-Point (FastPFP), 
which is capable of dealing with large graphs of over 
1,000 nodes. By using a new projected fixed-point 
method and a new partial doubly stochastic projection, 
our algorithm has time complexity O(n^) per iteration 
and space complexity O(n^). In addition to its scalability, 
our algorithm is easy to implement, robust, and able to 
match undirected weighted attributed graphs of differ- 
ent sizes. We also proved the linear convergence of the 
new projected fixed-point algorithm, based on the theory 
of Convex Projection [37] [38]. To the best of our knowl- 
edge, our algorithm is the only iterative graph matching 
algorithm with linear convergence guarantee. We con- 
duct extensive experiments on benchmark datasets of 
large graphs. FastPFP demonstrated better performance 
compared to previous state-of-the-art algorithms in most 
cases (only except for the graph isomorphism case, in 
which Umeyama's method is the fastest algorithm). In 
particular, in a PC, FastPFP is able to match two graphs 
of 1,000 nodes within a few seconds. 

The rest of this paper is organized as follows: In 
Section 2, we review previous work on continuous 
relaxation based graph matching algorithms. A sim- 
plied analysis of performance of different algorithms 
is given. In Section 3, we present our formulation of 
the graph matching problem. In Section 4, we introduce 
our FastPFP algorithm including the derivation, the new 
projection method, and convergence analysis. In Section 
5, we show extensive experiments conducted on various 
benchmark datasets in comparison with previous state- 
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of-the-art fast algorithms in large graph matching. In 
Section 6, we discuss FastPFP's comparison to FastGA, 
parameter sensitivity and limitations of our algorithm. 
Finally, we give concluding remarks in Section 7. 

2 Previous Work 

In this section, we review and analyze continuous relax- 
ation based graph matching algorithms. The main idea 
of continuous relaxation approach is to relax the original 
discrete optimization problem into a continuous one and 
solve it by continuous optimization techniques. For a 
simplified analysis, let us consider the performance of 
different algorithms based on continuous relaxation in 
matching two weighted graphs of n nodes. Relaxation 
Labeling (RL) (HI, EO), lEl, EH is an iterative updating 
algorithm with the assumption that the likelihood of two 
nodes to be corresponded can be inferred from their 
neighboring nodes. The algorithm has time complexity 
O(n^) per iteration. Graduated Assignment (GA) [23] is 
another iterative algorithm based on continuous relax- 
ation. It enjoys a clear optimization objective function 
and is widely regarded as a state-of-the-art algorithm 
in the field. Several techniques, including continuous 
relaxation, linear approximation, softmax and Sinkhorn 
normalization, are used to optimize this objective func- 
tion. Compared to RL, GA also has time complexity 
O(n^) per iteration but is more accurate and robust 1231 . 
Although the Linear Programming (LP) [24] approach is 
also based on linear approximation as GA, it has time 
complexity O(n^). 

The recently proposed new graph matching 
algorithms include Projections Onto Convex Sets 
(POCS) |25l. Spectral Matching (SM) [26], Spectral 
Matching with Affine Constraint (SMAC) [27], Integer 
Projected Fixed Point (IPFP) [29], and Reweighted 
Random Walks Matching (RRWM) ||30l. They all have 
time complexity O(n^) (per iteration), since they all 
require to construct and compute on a compatibility 
matrix or an association graph, which is an x 
matrix. Probabilistic Graph Matching (PGM) [28] 
requires the construction of a marginalization matrix, 
which takes O(n^) operations as well. Note that, 
although it is claimed that the compatibility matrix 
can be very sparse so that sparse matrix techniques 
can be used for efficient storage and computation, this 
is not the case for either weighted graphs or densely 
connected unweighted graphs. Path Following (PATH) 
algorithm [31] has time complexity O(n^) per iteration. 
However, it calls the Hungarian method [32] many times 
for each iteration. Although the Hungarian method has 
time complexity O(n^), the constant factor inside the 
big-O notation is large. In a PC, it takes more than 500 
seconds to apply Hungarian method to a 1000 x 1000 
non-sparse random matrix (See Fig. Therefore, these 
methods do not have good scalability with respect to 
graph size. Although some of the above algorithms 
have high matching accuracy, their high computational 
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Fig. 1. Runtime of Hungarian method applied to random 
non-sparse matrices. Matrix size ranges from 100 x 100 to 
1500 X 1500. 



costs prohibit them from matching large graphs (e.g. 
graphs of 1,000 nodes) within a reasonable time. 

Distinct from the above algorithms, there are three 
fast graph matching algorithms, all having time com- 
plexity O(n^) (per iteration): Umeyama's methods [33], 
FastGA [23 1 and Projected Gradient |34|. Umeyama's 
method uses eigendecomposition and Hungarian meth- 
ods only once, both of which have time complexity 
O(n^). In 1311 , it was shown that Umeyama's method 
is more than an order of magnitude faster than PATH. 
Despite of its efficiency, Umeyama's method is unrobust, 
especially when the two graphs to be matched are far 
away from isomorphic (see Section |5|. Moreover, it can 
only match graphs of equal sizes, which largely limits 
its use in practice. GA in general has time complexity 
O(n^) per iteration. However, by using a simple trick 
(see Appendix), GA can be directly reduced to have 
time complexity O(n^) per iteration. We call the reduced 
GA as FastGA. Although IPFP can also be reduced 
to O(n^) per iteration using a similar trick, IPFP calls 
the Hungarian method for each iteration and it takes 
5 10 iterations for IPFP to converge |29|. Therefore, 
the reduced IPFP is still not suitable for large graph 
matching. Projected Gradient (PG) is a classic method 
to solve constrained optimization problems including 
graph matching. However, PG can get stuck in poor local 
optima easily in solving the graph matching problem, as 
we observed empirically (Section |5|. The characteristics 
of different graph matching algorithms based on contin- 
uous relaxation are summarized in Table 1. 

The proposed new projected fixed-point algorithm for 
graph matching is substantially different from the pre- 
vious projected fixed-point method IPFP in three ways: 
(1) The projection of IPFP is onto a discrete domain, the 
space of partial permutation matrices, while the projec- 
tion of FastPFP is onto a continuous domain, the space 
of partial doubly stochastic matrices, by using a new 
projection method |40 | (Section [43) . (2) IPFP requires the 
expensive Hungarian method for each iteration, making 
it unsuitable for large graph matching. On the other 
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TABLE 1 

Characteristics of graph matching algorithms based on continuous relaxation 



Algorithm 


Time 
Complexity 


Space 
Complexity 


Able to Match Graphs 
of Different Sizes? 


Calling Hungarian method 
for each iteration? 


Convergence 
Guarantee 


Convergence 
Rate 


LP 


o(V) 


o(V) 


No 


No 


N/A 


N/A 


RL 


O(n^) /iteration 


0(n2) 


Yes 


No 


Yes 


Unknown 


GA, POCS 


O(n^) /iteration 


0{n^) 


Yes 


No 


Yes 


Unknown 


RRWM 


O(n^) /iteration 


0{n^) 


Yes 


No 


No 


N/A 


IPFP 


O(n^) /iteration 


0{n^) 


Yes 


Yes 


Yes 


Unknown 


SM, SMAC 




0(n4) 


Yes 


No 


N/A 


N/A 


PGM 




0(n2) 


Yes 


No 


N/A 


N/A 


PATH 


0(n'^) /iteration 


0(n2) 


Yes 


Yes 


Yes 


Unknown 


Umeyama 


0(n3) 


0(n2) 


No 


No 


N/A 


N/A 


FastGA, PG 


0(n'^) /iteration 


0(n2) 


Yes 


No 


Yes 


Unknown 


FastPFP 


0(n'^) /iteration 


0(n2) 


Yes 


No 


Yes 


Linear 



hand, each iteration of FastPFP only requires simple 
arithmetic operations (matrix addition and multiplica- 
tion). (3) IPFP is proved to converge at finite steps | 29| 
but the convergence rate is unknown. On the other hand, 
FastPFP is proved to converge at a linear rate. 

3 Problem formulation 

For two undirected graphs of size n and n' , denote 
their adjacency matrices (symmetric) by A and A' (bi- 
nary valued for unweighted graphs and real valued for 
weighted graphs) and attribute matrices by B and B' , 
respectively. Each row of B and B' is a /c-dimensional 
vector representing the attributes of a node. The size of 
A, A' , B and B' are n x n, n' x n' , n x k and n' x k, 
respectively. The objective function of graph matching is 

min - XA'X^Wl ^ \\\B - XB'\\l, (1) 

s.t. XI < = i,x e {o,ir><^', (2) 

where ||.||f is Frobenius matrix norm, A is a control 
parameter, and 1 is a vector with all its elements equal 
to one. In the left term can be interpreted as dissim- 
ilarity between edges and the right term as dissimilarity 
between nodes. Constraints in (2) enforces that X is 
a partial permutation matrix. The constant ^ is for 
convenience, to be seen later. We assume without losing 
generality that n > n' in the paper. The minimization 
problem ([l|(|2| is equivalent to 

max -tr{X^AXA)^Xtr{X^K), (3) 
JC z 

s.t. XI < i,x^i = i,x e {o,i}^><^', (4) 

where tr denotes the matrix trace and K denotes the 
nx n' matrix BB'^ (See Appendix for derivation). This 
problem is a Quadratic Assignment Problem (QAP) Il36l 
and is in general NP-hard. 



4 Algorithm 

Our proposed algorithm, FastPFP relies on continuous 
relaxation of the original discrete problem. The domain 
of the original discrete problem, the space of partial 
permutation matrices 

{X I XI < 1, X^l = 1, X G {0, 

is relaxed to onto the space of partial doubly stochastic 
matrices 

{X I XI < 1,X^1 = 1,X > 0}. 

A projected fixed-point method is then defined on the 
relaxed space. In the following subsections, we will give 
the derivation of our algorithm. 

4.1 Fast Projected Fixed-Point 

By relaxing the domain of the original QAP problem 
(|3|(|4| onto the space of partial doubly stochastic matri- 
ces, we have 

max \tr{X^AXA')^\tr{X^K), (5) 

S.t. XI < 1, X^l = 1, X > 0, (6) 

where X > denotes that all elements of X are nonneg- 
ative. To derive the algorithm, note that for 

/(X) = ]^tr{X^AXA') + \tr{X^K), (7) 

where X G the gradient of /(X) is 

V/(X) = AX A' + \K. (8) 

In this paper, we introduce a new projected fixed-point 
algorithm defined as 

X^^+i) = (1 - a)X^^) + aP,(V/(X(^))), (9) 

s.t. aG[0,l], (10) 

where Pd{-) is a partial doubly stochastic projection, 
defined as 

P^(X) = argmin||X-d||^, (11) 
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s.t. dl<l,(fl = l,d>0, (12) 

and a is the step size parameter. For a = 1, the algorithm 
is a straightforward projected fixed-point method. For 
< a < 1, the algorithm is a proportionally updated 
projected fixed-point method. The proportionality intro- 
duces more smoothness of the updating process, which 
could help stabilizing the algorithm (Section [6!2| ). We call 
this new algorithm Fast Projected Fixed-Point (FastPFP) 
due to its linear convergence guarantee (Section |4.4[ ). 
It also enjoys a property: if the initial state X^^'> is a 
partial doubly stochastic matrix, then each X^^"^ stays in 
the space of partial doubly stochastic matrices. This is 
because for two partial doubly stochastic matrices X and 
F, their convex combination 

Z = {l-a)X^aY, s.t. aG[0,l], (13) 

is another partial doubly stochastic matrix. 

4.2 Interpretation and related work 

In (|9|, let a = 1 and replace Pd{-) with the projection 
onto the space of partial permutation matrices, defined 
as 

PpermiX) = arg mm ||X - P\\%, (14) 

s.t. PI < 1, P^l = 1, P G {0, 1}^><^' , (15) 

which can be solved by the Hungarian method. Then (|9| 
becomes essentially a simplified IPFP 

= Ppe™(V/(X(*))). (16) 

The objective function (|3| is increased for each step 
of IPFP, as proved in |29|. Therefore, FastPFP can be 
interpreted as a softened version of IPFP. 

On the other hand, the Projected Gradient (PG) algo- 
rithm is defined as 

^(^+1) = Pa{X^'^ + aV/(X('))). (17) 

Although FastPFP and PG are algorithmically similar, in 
the graph matching problem, PG can get stuck in poor 
local optima easily. Moreover, PG usually runs signifi- 
cant slower than FastPFP. These experimental findings 
will be shown in detail in Section |5] And despite of 
the convergence guarantee of PG [34J, the convergence 
rate of the PG algorithm for solving the graph matching 
problem is unknown. 

In the next subsection, we will discuss how to obtain 
the solution of the projection Pd{-)- 

4.3 Partial Doubly Stochastic Projection 

The partial doubly stochastic projection can be converted 
to doubly stochastic projection by introducing slack vari- 
ables. Let Y be the slacked matrix of X 





fXu 


X12 


Xin' 


^l(n' + l) • 




Y = 


X21 


X22 


' X2n' 


^2(n' + l) • 


■ Yin 






Xn2 • 




^n(n' + l) 


^nn j 



(18) 



We write the above definition of Y compactly as 
^i:n,i:n' = X, where Yx,n,\:n' deuotes the matrix formed 
by the first n rows and the first n' columns of Y . 
To project a real nonnegative matrix onto the space of 
doubly stochastic matrices, the Sinkhorn method |39 | is 
usually used [231, |28l, ||30|, which normalizes each row 
and column of a matrix alternatively. The objective of 
the normalization is to find a nearest doubly stochastic 
matrix 1) to a matrix Y under the relative entropy 
measure. Zass and Shashua proposed another doubly 
stochastic projection |40|, which is originally to find 
a symmetric nearest doubly stochastic matrix D io Y 
under the Frobenius norm. However their method is also 
applicable to asymmetric doubly stochastic matrices, if 
Y is asymmetric. The projection of a matrix Y onto the 
space of doubly stochastic matrices under the Frobenius 
norm is the solution to the following problem 

Pz,(y) = argmm||y'-I)||^, (19) 

s.t. Dl = 1, D^l = 1,D>0. (20) 

With this doubly stochastic projection, we are able to 
prove the linear convergence of FastPFP in the next 
subsection. 

The doubly stochastic projection problem can be 
solved by successive projection. Define two sub- 
problems (projections) of ( p^pO) as 

Pi{Y) = argmin ||r - D\\f, s.t. Dl = 1, D^l = 1, 

(21) 

P2(F) = argmin Dili., s.t. D>0. (22) 
Both Pi and P2 have a closed-form solution 

1 1 ^ Y'l 1 1 

P,{Y) = ¥ + {-!+ - -Y)ll^ - -ll^Y, (23) 

Y + in 

P2{Y) = ^^. (24) 

The derivation is left in Appendix. The successive pro- 
jection works as follows: Pd{Y) = ...P2PiP2PiP2Pi{Y). 
The von Neumann successive projection lemma 1411 
guarantees the successive projection converges to Pd{Y). 
The solution to the partial doubly stochastic projec- 
tion under Fronbenius norm is therefore Pd{X) = 

PD{Y)l:n,l:n' • 

4.4 Convergence Analysis 

Now we state the convergence theorem of FastPFP. The 
proof of the theorem relies on Lemma 1. 

Lemma 1: Define Q = {X\X1 < = 1,X > 0}. 

For two real matrices Xi and X2, \\Pd{Xi) - Pd{X2)\\F < 
\\Xi-X2\\f- 

Proof: Since is a closed convex set, Pd(-) is a 
nonexpansive projection, which means ||Xi — X2||f ^ 
\\Pd{Xi) — Pd{X2)\\F' See [38 J for the proof of nonexpan- 
sivity of projections onto convex sets. □ 
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Theorem 1: Given real matrix X'^^\ a e [0, 1] and \\A<^ 
A' II 2 < e, the series 

= (1 - a)X^'^ + aPd{AX^'^A' + XK), (25) 

converges at rate 1 — a + ae. 

Proof: Denote vec{X), vec{Y), vec{K) and A^A^ by 
X, y, k and A, respectively, where vec{-) is the vectoriza- 
tion of a matrix and is the Kronecker product. Then 
Y(t) = AX^^^A' + \K is equivalent to y^^^ = Ax^^) + Ak. 
Thus 



- X(^) 11^ < (1 - a) • ||X(^) - X(^-i) 



(26) 



a-\\Pd{AX^'^A' ^\K) (27) 
P^(AX(^-i)A' + Ai^)||F. (28) 



Also 



\\Pd{AX^'^A' + Ai^) - Pd^AX^"^-^^ A' + Ai^)||F (29) 

< WAX^'^A' - AX^'-^^A'Wf = ||Ax(^) - Ax(^-^)||2 (30) 

<||A||2-||x(^)-x(^-i)||2 (31) 

<e-||X(^)-X(^-i)||F. (32) 

Thus 

-X(^)||7. 

□ 

Remarks. In practice, we do not have to scale A and A' 
such that II A (g) A'||2 < 1. We can let X = X/max(X) at 
the end of each iteration to prevent numerical instability. 

4.5 Discretization 

After the convergence of the projected fixed-point 
method, the resulting matrix X is discretized to 
obtain the partial permutation matrix P. We use the 
same discretization method in [26J, which is a greedy 
algorithm, instead of the expensive Hungarian method. 
This greedy discretization algorithms works as follows: 

Step 1. Initialize an n x zero-valued matrix P and a 
set L containing all its index 

Step 2. Given a matrix X, find the index (i*, j*) from 
L such that Xi^j^ = argmax(i .,)^^ X^^. Set Pi^j^ = 1. 
Remove all indices (i, j) in L that i = or j = j*. 
Step 3. Repeat step 2 until L is empty. Return P. 

The whole algorithm is summarized in Algorithm 1. 
Note that n > n' is assumed. The matching accuracy 
of FastPFP is sensitive to initialization. We recommend 
the initialization X = /nn', F = Onxn/ as we used 

in all our experiments. Step 9 is needed for numerical 
stability. For n = n' , regardless of fast and sparse ma- 
trix computation. Step 3 in Algorithm 1 requires O(n^) 
operations per iteration. Step 4 to Step 7 requires O(n^) 
operations per iteration. Step 8 requires O(n^) operations 
per iteration. Thus, the algorithm has time complexity 
O(n^) per iteration and space complexity O(n^). 



Algorithm 1: Fast Projected Fixed-Point 

Input: A,A',K 

Output: P 
1 Initialize X and Y 
1 repeat 



AX A' + \K 



1^ Yl 



Yl:n,l:n' 

repeat 

until Y converges; 
X = (l-a)X + aYi:n,i: 
X = X/max(X) 

10 until X converges', 

11 Discretize X to obtain P. 



I - ^Y)H^ - ^ll^Y 

n / n 



5 Experiments 

In this section, we evaluate the performance in large 
graph matching of four fast algorithms: FastPFP, Pro- 
jected Gradient (PG) jTzj, FastGA and Umeyama's 
method. All of them have time complexity O(n^) (per 
iteration). Some other graph matching algorithms were 
not compared due to their inability to match large graphs 
(e.g. graphs of 1000 nodes) in a reasonable time. 

For FastGA, FastPFP, and PG, the same greedy dis- 
cretization described in the above section was used. 
The convergence criteria are: max(|X'^^+^^ — X'^^^l) < ei 
or more than /q iterations for the first loop (projected 
fixed-point updating for FastPFP or soft-assignment for 
FastGA) and max 

_ < e2 or more than h 

iterations for the second loop (doubly stochastic projec- 
tion or normalization). For the fairness of comparison, 
we tuned and fixed the parameters of each algorithm 
throughout all our experiments. The parameters were 
well-tuned to achieve the best performance taking both 
accuracy and efficiency into consideration. For FastPFP, 
we fixed a = 0.5 in all experiments. The sensitivity study 
of this parameter will be given in the next section. Other 
parameters were also tried to validate the experimental 
results. All our experiments were implemented in MAT- 
LAB in a 3 GHz Intel Core2 PC. 

Remarks. (1) In some datasets, the experiments of 
Umeyama's method were not conducted due to its in- 
ability of matching graphs of different sizes. (2) In the 
weighted graphs matching experiments, graphs are fully 
connected. Each node represents a keypoint and each 
edge represents Euclidean distance between two key- 
points. We tried Delaunay triangularization to sparsify 
the graphs. But the matching quality was poor after 
applying the sparsification. Therefore, no sparsification 
was used in the experiments, for simplicity and better 
matching results. 

5.1 Synthetic Random Graphs 

In this set of experiments, unweighted random graphs 
were generated uniformly with 50% connectivity. The 
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(a) No noise (graph isomorphism) 



(b) n edges edition 
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(c) 10% nodes deletion (subgraph isomorphism) 

Fig. 2. Synthetic random graph matching 



(d) 10% nodes deletion + n edges edition 



size of the graphs ranges from 100 to 1,500. Four different 
types of graph matching were tested: two isomorphic 
graphs with (a) no noise (graph isomorphism), (b) n 
edges edition (edge edition means if Aij = 1, set 
Aij = and vice versa), (c) 10% nodes deletion (sub- 
graph isomorphism), (d) n edges edition + 10% nodes 
deletion. The runtime and matching error (measure by 
\\A-XA^X^\\l - \\A-XgtA^X^^\\1, where Xgt is the 
ground truth matching matrix) were recorded, as shown 
in Fig. H 



As observed from Fig. |2| FastPFP and FastGA were 
able to find the global optimal solutions in all the above 
graph matching tasks. But FastPFP is about 3^5 
times faster than FastGA. PG were stuck in poor local 
minimum in most cases. Umeyama's method achieved 
better performance than other algorithms for the iso- 
morphic graph matching in terms of both speed and 
accuracy. However, its performance was poor in match- 
ing non-isomorphic graphs. Roughly speaking, the more 
isomorphic of two graphs, the better matching results 
Umeyama's method had. In the left part of Fig.|2|b) , the 
matching error increases with the number of nodes from 
100 to 500. This is natural since the numeric range of the 
matching error increases with the number of nodes given 
the same percent of connectivity. From 600 nodes on, 
the matching error decreases with the number of nodes. 
This is because the graphs are closer to isomorphic as the 
number of nodes increases since there are only n edges 
edition for total 0{n^) edges. 



5.2 Graphs From Image Sequences 

In this set of experiments, weighted graphs were con- 
structed from the CMU House sequence image^ which 
has been wildly used for testing graph matching al- 
gorithms. In the previous literature, only small graphs 
(30 nodes) were extracted from the images for match- 
ing. To demonstrate scalability of FastPFP, we applied 
SIFT feature descriptor [42] to all CMU House sequence 
images. Roughly 800 keypoints were extracted for each 
image and the number of keypoints varies from image 
to image. All keypoints in an image were used for 
the construction of a graph. For each image, a graph 
whose nodes represent keypoints and edges (matrix ^4) 
represent Euclidean distance between keypoints in that 
image was constructed. For a visualization, we show in 
Fig. |3] feature correspondence for image vs. image 110. 
Results of PG, FastGA and FastPFP were shown. 

We matched the first image and subsequent images, 
i.e. image vs. image 1, image vs. image 2, and 
image vs. image 110. The runtime and matching error 
(measured as HA-XA^X^Hl.) were recorded for different 
algorithms. The results are shown in Fig. |4] FastPFP is 
about 6 times faster than FastGA, with even much lower 
matching error. 

5.3 Graphs From Real Images 

In this set of experiments, attributed and weighted 
graphs were constructed from real images, chosen from 
the four image sets: Graffiti, Wall, Bark and Boal|^ The 
graphs were constructed as follows. We first applied 
SIFT descriptor to extract 128-dimensional features for 



http:/ / vasc.ri.cmu.edu/idb/html/ motion/ 
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i/affine/ 



(a) SIFT features (820 nodes vs. 772 nodes) 




(b) From top to bottom: PG, FastGA and FastPFP 




(c) SIFT features (820 nodes vs. 772 nodes) 




(d) From top to bottom: PG, FastGA and FastPFP 

Fig. 3. Graphs from CMU house sequence matching. For (c)(f), Left column: 10% matching displayed. Middle column 
20% matching displayed. Right column: 100% matching displayed. 
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Fig. 4. CMU house sequence 



TABLE 2 

Graphs from real images: runtime (sec) 



Graph pair 


PG 


FastGA 


FastPFP 


Graffiti 


51.60 


87.89 


13.52 


Wall 


19.94 


32.25 


4.92 


Bark 


2.84 


2.72 


0.52 


Boat 


5.59 


9.06 


1.36 



TABLE 3 

Graphs from real images: matching error 



Graph pair 


PG 


FastGA 


FastPFP 


Graffiti 


2.75 xlO"^ 


3.98 xlO"^ 


5.92 X 10^ 


Wall 


1.31 xlO^ 


1.01 xlO^ 


2.39 X 10^ 


Bark 


3.07x10^ 


2.21 X 10^ 


2.21 X 10^ 


Boat 


4.77x10^ 


3.18x10^ 


1.80 X 10^ 



each image. Then a simple feature selection was per- 
formed as follows: a feature of an image was selected 
if it has the similarity (inner product) to all features 
of its counterpart image above a threshold. From the 
selected features, a graph was constructed. Each node 
represents a keypoint, each node attribute represents 
a SIFT 128-dimensional feature vector, and each edge 
represents the Euclidean distance between nodes. Each 
element Kij of K is set to be the inner product of 
two SIFT feature vectors i and j. See Fig. |6] for visu- 
alization of the matching results for PG, FastGA and 
FastPFP. The experimental results are listed in Table |2] 
and |3] The runtime and matching error (measured as 
I \\A - XA'X^Wl + A||5 - XB'Wl) were recorded. A = 1 
is used. The sensitivity study of parameter A is given 
in Section |6.2[ Again, FastPFP demonstrated the best 
performance: it is about 3 ~ 6 times faster than FastGA 
while usually achieving a much lower matching error. 



5.4 Graphs From 3D Points 

In this set of experments, weighted graphs were con- 
structed from 3D points. The Face 50 and Stanford 



TABLE 4 

Graphs from 3D points: runtime (sec) 



Graph pair 


PG 


FastGA 


Umeyama 


FastPFP 


Face pair 1 
Face pair 2 
Stanford Bunny 


2.90 
3.06 
21.74 


15.05 
15.53 
118.55 


18.40 
18.23 
1027.02 


0.69 
0.69 
8.75 




TABLE 5 






Graphs from 3D points: matching error (10^) 


Graph pair 


PG 


FastGA 


Umeyama 


FastPFP 


Face pair 1 
Face pair 2 
Stanford Bunny 


6.30 
6.23 
29.1 


14.3 
16.2 
37.4 


12.5 
10.6 
32.7 


0.78 
1.64 
9.55 



Bunny 3D points dataset^ were used. In the Face 50 
dataset, each face has 392 points, all of which were used 
for the construction of graphs. In the Stanford Bunny 
3D data, each set of 3D points was downsampled to 
1022 points. Graphs, whose edge represents Euclidean 
distance between two points, were constructed. The 
visualization of matching results of PG, FastGA, FastPFP 
and Umeyama's method are shown in Fig. U\ and |8] 
The matching error (measured as \\A - XA^X^.) and 
runtime are shown in Table [4] and [5] Again, FastPFP 
outperformed the others in both speed and accuracy. 
FastPFP is about 15 ^ 20 times faster than FastGA. 

6 Discussion 

6.1 Comparison To FastGA 

Although FastGA and FastPFP both have time com- 
plexity O(n^) per iteration, FastPFP is empirically about 
3 ~ 20 times faster than FastGA. This is mainly due to 
the slow convergence of FastGA. Despite the theoretical 
work on convergence of FastGA [35], empirically FastGA 
does not have good convergence speed, as observed in 
the experiments. To see that, we show the magnitude 

3. Both datasets can be downloaded in https:/ / sites.google.com/ j 
site / myronenko / 1 



(a) SIFT features (1438 nodes vs. 1443 nodes) 




(b) From top to bottom: PG, FastGA and FastPFP 



(c) SIFT features (961 nodes vs. 960 nodes) 




(d) From top to bottom: PG, FastGA and FastPFP 

Fig. 5. Graphs from real images matching. For (c)(f), left column: 10% matching displayed. Middle column: 20% 
matching displayed. Right column: 100% matching displayed. 



(a) SIFT features (393 nodes vs. 391 nodes) 




(b) From top to bottom: PG, FastGA and FastPFP 




(c) SIFT features (542 nodes vs. 537 nodes) 




(d) From top to bottom: PG, FastGA and FastPFP 



Fig. 6. Graphs from real images matching. For (c)(f), left column: 10% matching displayed. Middle column: 20% 
matching displayed. Right column: 100% matching displayed. 




(b) PG 



(c) FastGA 




(d) Umeyama 




(e) FastPFP (Ours) 



Fig. 7. Graph from 3D points matching 
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(c) FastGA 





(d) Umeyama 



(e) FastPFP (Ours) 

Fig. 8. Graph from 3D points matching: Stanford Bunny. For (b)-(e), left column: 10% matching displayed. Middle 
column: 20% matching displayed. Right column: 100% matching displayed. 
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5 10 15 
Number of iterations 

(a) image sequence sequence 



5 10 15 
Number of iterations 

(b) Real image 



Fig. 9. Convergence of FastPFP and FastGA. The vertical 
axis represents the value of ei = max(|X'^^+^) - X^^^\). 
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(b) Real image 

Fig. 1 0. Performance for different values of a 



of parameter ei = max(|X*^^+^^ — X*^^^|) for image se- 
quence (image vs. image 110) and real image (Graffiti) 
for 20 iterations in Fig. [9] Besides good convergence, 
FastPFP only requires simple arithmetic operations (ad- 
dition, multiplication and maximum element search) 
while FastGA requires additionally exponentiation and 
logarithm (for numerically stable implementation). 

6.2 Parameter Sensitivity 

The performance of FastPFP is not very sensitive to the 
step size parameter a, for < a < 1. In Fig. 10 the 
runtime and matching error of FastPFP in the image 
sequence (image vs. image 110) and the real image 
experiment (Graffiti) are shown for different values of 
a, ranging from 0.1 to 1. The performance of FastPFP is 
unstable for a = 1. 

For attributed graph matching, the matching quality 
is not very sensitive to the parameter A. We show the 
matching result for A = oo (which means the matching 
objective function becomes \\B — XB'Wj^ and is solved 



(a) A = oo 




(b) A = 



(C) A: 




(d) A = 10 

Fig. 1 1 . Performance for different values of A. 1 0% of the 
matching are displayed. 



by Hungarian method), A = and A = 10. The results 
are similar except for A = oo. The matching results are 
shown in Fig. 
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6.3 Limitations 

Despite the encouraging experimental results, our algo- 
rithm is derived heuristically, as many other approxi- 
mate graph matching algorithms such as GA (FastGA), 
Umeyama's method, SM and RWRW. For these algo- 
rithms including ours, the optimality of the solutions is 
not theoretically guaranteed. Also, our algorithm does 
not use an edge compatibility matrix, due to its high 
computational cost. Consequently, another limitation of 
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our algorithm is its inability to handle arbitrary edge 
relationship (some graphs have vector edge attributes), 
unlike those algorithms based an edge compatibility 
matrix. 

7 Conclusion 

We proposed a new fast graph matching algorithm 
based on a new projected fixed-point method, suitable 
for large graph matching. Extensive experiments were 
conducted to demonstrate the strength of our algorithm 
over previous state-of-the-art algorithms. Future work 
includes theoretical analysis, extension and applications 
of our algorithm. 
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Appendix 

FastGA 

In general, GA also has time complexity O(n^) per itera- 
tion. However, when the compatibility matrix C = A^A' 
(which was mentioned in the original paper of GA |23|) 
is used, GA has only time complexity O(n^) per iteration. 
To see this, consider the most time-consuming step of 
GA, which has time complexity O(n^) 

n n 

b=i j=i 

(34) 

where X is the soft-assignment matrix. If Caibj = ^ij^'ah' 
then can be compactly written as 



(35) 



which has time complexity O(n^). 
Problem Formulation 

By expanding the Fronbenius norm, we have 

\\\A-XA'X^\\l^\\\B-XB'\\l 



Successive Projection 

The derivation of (23|(24| is due to fiOl, with slight mod- 
ification in here. We first consider Pi . The Lagrangian of 
(|21) is 

L{D, ui,U2) = tr{D^D - 2X^D) (43) 
- uJ{Dl - 1) - u^{D^l - 1). (44) 

Let ui = U2 = u and set the derivative with respect to D 
to zero, we have 



D = X^ur 



Multiplies by 1 on both sides: u = {nl + 11 



(45) 

1 = (l/nMI - 



Combining the fact that (nl + 11^) 
(l/2n)ll^), we obtain The derivation of (|24|) is 

straightforward and omitted. 
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(36) 



^tr{AA^) + ^tr{X^XAX^XA^) - tr{X^AXA') 
^ ^ (37) 



-^\tr{BB^) + \tr{X^XB'B'^) - 2\tr{X^BB'), (38) 

due to the invariance of the matrix trace under cyclic 
permutation. Since X is a partial permutation matrix and 
n > n', X^X = 1. And the first two terms of ^ and 
([38) are constants. Therefore, the minimization problem 
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min -P - XA'X^Wl + X\\B - XB'\\l, (39) 



s.t. XI < 1, X^l = 1, X G {0, , (40) 

is equivalent to 

max -tr{X^AXA)^Xtr{X^K), (41) 
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s.t. XI < 1, X^l = 1, X G {0, . (42) 



